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T^lj- ■ A new method for solving the configuration-space Faddeev equations for 

. elastic p-d scattering below the deuteron-breakup threshold is described. Nu- 

merical solutions that demonstrate the convergence and accuracy of the method 
are given. The number of channels and the value of the matching radius re- 
quired to obtain an accurate solution are also investigated. These calculations 
^ , demonstrate that this method can efficiently solve the large matrix equations 

^ ' required for the three-body scattering problem. 

c ■ 



Abstract. 



1 Introduction 

The three-nucleon system is a useful tool for investigating our understanding of 
the nuclear force. Since it is possible to obtain accurate numerical solutions for 
the three-body system, detailed comparisons of the model Hamiltonian with the 
experimental data can be performed. These comparisons provide stringent tests 
of the two-body interactions obtained by fitting the two-body data, and give a 
means of studying the three-nucleon force. Accurate bound-state calculations ||l| 
have shown that realistic two-body interactions yield binding energies for the 
triton and ^He that are less than the experimental values, and one must add 
a three-body interaction to obtain the correct value. The bound-state system 
unfortunately provides limited information about the three-body force. Most of 
the bound-state properties scale with the binding energy [^. That is, different 
interactions that give the same trinucleon binding energy predict nearly the same 
charge radii, asymptotic normalization constants, magnetic moments, etc. 

Initially it was believed that the nucleon-deuteron scattering problem offered 
an opportunity to explore our understanding of the three-nucleon force; however, 
it has been shown Q that most of the scattering data can be reproduced at a good 
level with only two-body interactions. For low-energy scattering the effects of the 
three-nucleon force are usually small; there are nevertheless some discrepancies 
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(such as the nucleon analyzing power Ay and the deuteron analyzing power iTn) 
that could be sensitive to these interactions. Since much of the experimental data 
is for p-d scattering at low energies, where the Coulomb interaction cannot be 
neglected, it is necessary to solve the scattering equations for this case as well 
as for the n-d case. 

The Pisa group has shown that it is possible to obtain accurate solutions 
for the p-d scattering equations using the Pair-correlated Hyperspherical Har- 
monic (PHH) basis to expand the wave function, and solving for the S-matrix by 
using the complex form of the Kohn variational principle. In this paper we de- 
scribe an alternate method for solving the configuration-space equations. While 
the method is similar to the one used in our previous calculations of the scatter- 
ing length Q, a new iterative technique for solving the large matrix equations 
is presented. In addition to demonstrating the convergence and accuracy of this 
technique, we provide a detailed study of the number of angular-momentum and 
spin-isospin states required to obtain an accurate solution. Since the question 
of how large a matching radius is required in configuration space to obtain an 
accurate solution has been a concern about this method, we also investigate the 
convergence of the method as a function of this parameter. While the present 
work is limited to scattering below the three-body breakup threshold, the method 
can be extended to higher energies using the boundary conditions discussed in 
Ref. §. 

The new numerical method is described in the next section, and the conver- 
gence of the procedure is illustrated in Sect. 3. 



2 Numerical Methods 

We use the Jacobi coordinates 

Xi = - rk (2.1) 

and 

yi = ^i^j + n)-^i, (2-2) 

where i, j, and k imply cyclic permutation. The Hamiltonian in the center-of- 
mass frame is 

if = r + F(xi) + F(x2) + F(x3) + yc, (2.3) 

where 

[l+r,(j)][l+r,(fc)] 
Vc = l^- (2.4) 

1=1 ' 

is the sum over the two-body Coulomb interactions between the pairs and V{xi) 
is the nucleon-nucleon potential between the pair j and k. For an incident nu- 
cleon with kinetic energy E'lab in the laboratory frame, the total energy in the 
center-of-mass frame is the sum of the deuteron binding energy and the kinetic 
energy of the incident particle in the center-of-mass frame, Ecm = f-Eiab. The 
corresponding wave numbers are given by 

E-E, + E,^--^ + ^ = -^, (2.5) 
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where M is the nucleon mass, k is the deuteron bound-state wave number, and 
q is the wave number of the nucleon or deuteron in the center-of-mass frame. 
Writing the total wave function ^ as the sum of the three Faddeev amplitudes 

^ = *i(xi,yi) + *2(x2,y2) + *3(x3,y3) , (2.6) 

the Schrodinger equation can be decomposed into the three Faddeev equations 

[T + y(x,) + Vc- E]^i = -V{^i){^j + *fc) . (2.7) 

Adding these three equations gives the Schrodinger equation. For three identical 
particles the all have the same functional form, and we need to solve only the 
i = 1 equation. To solve this equation, we make a partial-wave expansion of the 
Faddeev amplitudes using the j — J coupling scheme. We write 

*i(xi,yi) = 2^ , (2-8) 

Q, Xiyi I 

where 

|a). = I [{la, Sa)ja, {La, Sa)Ja] JM; {ta,T^)TMT) (2.9) 

is the orbital angular momentum and spin-isospin state function for the different 
channels, and i indicates that the order of coupling is {{j,k)^i\ for cyclic values 
of i, j, and k. The relative orbital angular momentum of the j-k pair is la, the 
spin of the pair is s^, and the total angular momentum of the pair is ja. The 
orbital angular momentum of particle i relative to the center-of-mass of the j-k 
pair is La, Sa is the spin of particle i, and J a is the total angular momentum 
of the particle. Finally, J is the total angular momentum of the system. For the 
isospin function, ta is the isospin of the j-k pair, is the isospin of particle z, 
and T is total isospin of the system. For nucleons Sa and Ta are both 1/2. 

The Faddeev equation for »I'i(xi,yi) can be reduced to a set of coupled 
partial differential equations for the reduced amplitudes V'a by taking the inner 
product with the state functions la)^. After multiplying the equations by M/^^ 
and transforming to the hyperspherical coordinates defined by 

x = pcos9 (2.10) 

and 

y = ^psme, (2.11) 

the resulting equations arc 

(A, - /C2) ^a{p, e) - E (^"«' + V^a') ^a'{p, 0) 

a' 

- E E /' ^a'a"{0, e')i^a"{p, 0') dO' = 0, (2.12) 



a' a" 



where 



A _ 5' 1 5 1 la{la + 1) La{La + 1) 



dp^ p dp p^ 86"^ p^ cos^ 9 p^ sin 
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Vaa' = ^ia\Vi^i)\a') (2.13) 



and 



E /_ K^'a"{d,o')^Pc,^>{p,9')de' = xmy\^2{^2,y2) + M^3,y3))- (2.i4) 

„ // JO 

a 

For three equal-mass particles the limits for the integration with respect to the 
hyperspherical angle 9' are given by 



TT 



and 

TT 



6+ = ^ 
2 
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where < (9 < 7r/2. 

For the numerical solution of the coupled partial differential equations, we 
choose a value, Pmax, for the maximum value of the p variable, and solve the 
equations for < p < Pmax- The equations have a unique solution only if one 
specifies the boundary conditions on the closed boundary of this region of the 
p — 9 plane. Since the are reduced wave functions they must vanish as x or y 
goes to zero; thus, we have the boundary conditions 

V'a(0, 9) = ^Po^ip, 0) = Mp, ^/2) = (2.15) 

for all channels. The boundary conditions for pmax depend on whether the chan- 
nel is open or closed. For energies below the three-body-breakup threshold, the 
channels that correspond to elastic scattering are the only open ones. The other 
channels, which are virtual breakup channels, must decrease exponentially for 
large values of p. The channels whose values for /„, s^, and ja correspond to the 
deuteron values are open, and they have the asymptotic form 

'(l^a{xi,yi) — > 4>a{xi,yi) + ^a{xi,yi) , (2.16) 

where (i)a{xi,yi) is the component of the incident wave for the channel a and 
^a{xi,yi) is the corresponding component of the scattered wave. For p-d scat- 
tering 0a(a;i,yi) is given by 

(f>a{xi,yi) = FLSv,Qyi)ui^{xi) , (2.17) 

where FL^{ri,qyi) is the regular Coulomb wave function with rj = 2Me^ /3h'^q 
and ui^ is the la component of the reduced deuteron bound-state wave function. 
For scattering energies below the breakup threshold the K-matrix form of the 
scattering equations is real, and using this form of the equations simplifies the 
numerical calculations since the flaixi,yi) are real functions. For a deuteron 
with ja = 1 and a fixed value of the total angular momentum J and parity there 
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are at most three open channels, and the equations must be solved for incident 
waves in each of these channels. For an incident wave in channel a', the scattered 
wave in the K-matrix formulation has the form 



(2.18) 



where GL^jj^qyi) is the irregular Coulomb wave function. Given the K-matrix 
one can determine the S-matrix using 



l + iK 
l-iK 



(2.19) 



For the numerical calculations we write the reduced Faddeev amplitudes in 
the form 

) yi ) — ^aa' 4'ct {xi,yi) + xa{p,9), (2.20) 

for the case with an incident deuteron in the channels with a fixed value of Jq,/. 
We define the function Fa{p, 6) (which is not the regular Coulomb function used 
m Eq. ( p^ ) by writing 



Xa{xi,yi) = Fa{p,9)e 



-Kp 



for the closed channels, and 



Xa[Xl,yi) = Ui^ 



X\ 



(2.21) 



(2.22) 



for the open channels, where the deuteron wave function has been factored out 
of the elastic channels to make Foi{p,9) a smoother function. The xi in the 
denominator of Eq. ( 2.22 ) is included to force Fa{p, 6) to vanish when xi is zero. 
The boundary condition for closed channels is that Xaixi, yi) go to zero for large 
values of p, which can be implemented by the condition 



9F„ 



dp 



0: 



(2.23) 



P — Pmax 

that is, Fa{p,0) is a constant for large values of p. For the open channels 

Fa{xi,yi) n r ^ 

which is implemented by the boundary condition 
OF, 



(2.24) 



dp 



P — pmax 



cos 9 h — sm d— ; — 

xi 2 Gl^ dyi 



(2.25) 



To solve the equations for a given value of the total angular momentum, J, 
and parity, we truncate the number of channels by choosing a maximum value 
for the angular momentum, Jq, of the interacting pair and solve the differential 
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equation for an incident wave in each of the open channels. Substituting Eq. 
(1^) into Eq (|2l2|) gives 



a' 

-E^^^'ET Ko,'o."{0,e')xa"{p,e')d9' (2.26) 

a' a" 

= ^apY{v%' - Uj5i3pi)(t)p'{p,9) + Yvaa'Yl I ^a' f3' {0 , 9')(t)fj' {p, 9') dO' , 
P' a' 13' •'^ 

where the sum over (3' is over the channels that have an incident deuteron bound 
state in the asymptotic region, and u) = (M/^^)e^/?/i. Finally, substitution of 
Eqs. (|]2|) and (|23^ ) into (|236D gives a set of coupled partial differential equa- 



tions for the Fa{p, 9). To solve these equations we expand Fa{p, 9) in a complete 
set of basis functions that is the tensor product of Hermite splines for the p and 
9 variables. We write 

Np Ng 

Fa{p,9)=J2T.<^^(p>ji(^)^ (2-27) 
i=i j=i 

and use the orthogonal collocation method to determine the afj] that is, we 
require ( |2.27| ) to satisfy the differential equation at the collocation points {pi, 9m) 
for 1 = 1,..., Np and m = 1, . . . , Nq. This gives a matrix equation whose columns 
are labeled by the values of i, j, and a, and the rows are labeled by the values of Z, 
m, and a' . For an accurate solution to the equations the number of the expansion 
coefficients, a^-, can exceed several hundred thousand, and the resulting matrix 
equation is too large to invert directly. From Eq. ( 2.26| ) and the local nature of 



the splines, one can see that the matrix will have a block-diagonal structure. 
Since the tensor component of the nuclear force couples at most two channels, 
the matrix for the first two terms on the left-hand side of Eq. ( 2.26| ) will have a 



much smaller bandwidth than the third term, which couples all of the channels. 
Thus we write the matrix equation in the form 

{A-B)a = b, (2.28) 

where the matrix A can be inverted by standard methods, while the matrix B 
with a much larger bandwidth requires the use of an iterative procedure. 

Since standard iterative methods such as the Lanczos algorithm or Fade 
approximants required many iterations to converge to an accurate solution for 
the matrix equation, we derive a new algorithm to solve the matrix equation. 



Rewriting ( 2.28 ) in the form 

(l- BA-^^{Aa) = h, (2.29) 

we find an approximate solution for Aa using a set of orthogonal basis vectors 
constructed by the Gram-Schmidt procedure. Starting with uq = h/VWh, we 
iterate Eq. ( |2.29| ) to generate the basis vectors. Given a set of basis vectors, Ui 
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for i = 0,...,A^j — l,a new basis vector is constructed by first generating the 
vector 

VN, = (l - BA-^) UN,-i , (2.30) 

and then using the Gram-Schmidt procedure to generate a uat. which is nor- 
malized and orthogonal to the other basis vectors. After Ni iterations, the basis 
vectors Ui are used to find an approximate solution for Aa by writing 

Aa^ J2 ^i'^i ■ (2-31) 

i=0 

Substituting this approximation into Eq. ( p. 29 ) and using Eq. ( p.SOD , gives 

J2ciVi+i = b. (2.32) 

i=0 

Taking the inner product with Ui gives the matrix equation 

^f^i+i^i = ufb, (2.33) 

j=0 

which can be solved for the Cj. Multiplying the approximation for Aa by A~^ gives 
the approximate solution for a. Using this approximate solution, the Fa(p, 6) can 
be evaluated using Eq. ( 2.27 ), and the values of K^^a' can be determined using 
Eq. ( 2.24 ). The procedure is repeated until a solution of the desired accuracy is 
found. Normally 10-20 iterations are required to obtain an accurate solution. 



3 Numerical Results 



To demonstrate the convergence of the iterative method, we consider the case 
of p-d scattering for the complete AV18 two-body potential 0] with fi? /M = 
41.47108 MeVfm^. The K-matrix elements are for the J-j channel states defined 
in Eq. ( p. 91) , which are different than the channel states used for the n-d and p- 
d 1^ benchmark calculations. Our results for the K-matrices can be converted to 
the channel scheme using Eq. (2.26) in Ref. |^]. The incident state has T = 1/2, 
and we neglect the isospin mixing, which has been shown to be small for low- 
energy scattering [0]. Thus, for a given value of the total angular momentum, J, 
and parity, the possible channel states in Eq. ( p. 9]) are determined by the values 
of ja and Jet ■ For our calculations the number of channel states is determined by 
keeping all states with values of ja up to jmax- In addition, the equations must 
be solved numerically in the region < p < /Omax, where pmax must be large 
enough that the asymptotic boundary conditions are an accurate approximation 
to the scattered wave function. It is shown below that using jmax = 10 and 
Pmax = 90 fm gives accurate numerical solutions. To illustrate the convergence of 
the iterative solution to the numerical equations, we show in Table |^ the values 
for the K-matrix elements for the J = 1/2"*" scattering at E'lab = 1.0 MeV. For 
this case the indices 1 and 2 refer to the states with (Lq,, J^) equal to (0, 1/2) and 



8 



Low Energy Proton-Deuteron Scattering in Configuration Space 



Table 1. Convergence of the iterative solution for the 
K-matrix for J ~ p-d scattering at iJiab = 1.0 MeV 
evaluated using jmax = 10 and pmax = 90 fm. 



i 


Kii X 10^ 


Ki2 X lO'' 


K21 X 10^ 


K22 X 10^ 


3 


-1.3847 


-9.2557 


-2.5439 


-1.3052 


4 


-2.7632 


-2.5891 


-2.7734 


-1.2797 


5 


-2.4279 


-3.8322 


-3.5321 


-1.3369 


6 


-2.1698 


-4.3883 


-3.9463 


-1.3690 


7 


-2.3100 


-3.8242 


-4.0234 


-1.3701 


8 


-2.3003 


-4.0095 


-3.9951 


-1.3734 


9 


-2.3082 


-3.9661 


-3.9724 


-1.3734 


10 


-2.3092 


-3.9520 


-3.9597 


-1.3720 


11 


-2.3073 


-3.9685 


-3.9645 


-1.3721 


12 


-2.3076 


-3.9654 


-3.9653 


-1.3722 


13 


-2.3076 


-3.9641 


-3.9649 


-1.3722 


14 


-2.3077 


-3.9645 


-3.9650 


-1.3722 


15 


-2.3076 


-3.9645 


-3.9650 


-1.3722 


16 


-2.3076 


-3.9645 


-3.9650 


-1.3722 



(2, 3/2), respectively. The values for the Kij are calculated using the approximate 
solution obtained from Eq. ( 2.31 ). 

To show the accuracy of the solution using a truncated set of channel states, 
we list in Table |2| the values for the K-matrix elements calculated using different 
values for jmax for two different energies. The value of Nc is the actual number 
of channel states. Since the correct K-matrix must be symmetric, one test of the 
accuracy of the numerical solution is the difference between K^j and Kji. As the 
number of channels is increased, this difference becomes smaller. From Table |2| 
one can see that a value of jmax = 10 gives an accurate approximation for the 
K-matrix. 



Table 2. Values of the K-matrix for J ~ 2 P~d scattering evaluated using 
different numbers of basis states and pmax ~ 90 fm. 



1.0 





Jmax 


Kii X 10' 


K12 X 10" 


K21 X 10^ 


K22 X 10 


34 


4 


-2.3131 


-.38871 


-.39693 


-1.3737 


50 


6 


-2.3081 


-.39424 


-.39688 


-1.3727 


66 


8 


-2.3078 


-.39587 


-.39663 


-1.3724 


82 


10 


-2.3076 


-.39645 


-.39650 


-1.3722 



3.0 



34 


4 


-6.2352 


-1.0738 


-1.0859 


-6.2915 


50 


6 


-6.2261 


-1.0857 


-1.0872 


-6.2888 


66 


8 


-6.2253 


-1.0882 


-1.0870 


-6.2878 


82 


10 


-6.2250 


-1.0887 


-1.0870 


-6.2874 



Finally in Table H W6 show th.6 results for different va;lues of Pmax- 

A value of 

Pmax = 90 fm is more than enough to obtain an accurate solution to the scattering 
equations. 
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Table 3. Values of the K-matrix for J = ^ p-d scattering as a 
function of the laboratory energy in MeV evaluated using jmax ~ 10 
and different values of pmax in fm. 





pmax 


7^11 X 10^ 


Ki2 X 10^ 


K21 X 10^ 


A'22 X 10" 




30 


-2.2966 


-.37960 


-.39970 


-1.1755 


1.0 


50 


-2.3076 


-.39681 


-.39681 


-1.3762 














70 


-2.3076 


-.39652 


-.39657 


-1.3723 




90 


-2.3076 


-.39645 


-.39650 


-1.3722 




30 


-6.2440 


-1.0919 


-1.0972 


-6.5249 


3.0 


50 


-6.2212 


-1.0887 


-1.0860 


-6.2867 














70 


-6.2243 


-1.0888 


-1.0870 


-6.2878 




90 


-6.2250 


-1.0887 


-1.0870 


-6.2874 



In Tables ^ and ^ we show the convergence of the J = l/2~ state for a 
laboratory energy of 1.0 MeV. For this case the indices 1 and 2 refer to the 
states with (L^, Jq,) equal to (1, 1/2) and (1,3/2), respectively. The convergence 
is very similar to the J = 1/2+ case. 



Table 4. Values of the K-matrix for J ~ ^ p-d scattering at 
Slab ~ 1.0 MeV evaluated using different numbers of basis states 
and pmax = 90 fm. 





Jmax 


All X 10' 


K12 X 10^ 


K21 X 10^ 


K22 X 10^ 


34 


4 


1.4461 


-6.0489 


-6.0613 


-4.0558 


50 


6 


1.4451 


-6.0495 


-6.0542 


-4.0581 


66 


8 


1.4446 


-6.0497 


-6.0517 


-4.0589 


82 


10 


1.4445 


-6.0499 


-6.0508 


-4.0592 



Table 5. Values of the K-matrix for J ~ ^ p-d scattering 
at Blab = 1.0 MeV evaluated = 10 and different 

values of Pmax. 



Pmax 


All X 10' 


A12 X 10" 


A21 X 10^ 


A22 X 10" 


30 


1.4277 


-5.9366 


-5.9524 


-3.9590 


50 


1.4448 


-6.0514 


-6.0521 


-4.0616 


70 


1.4445 


-6.0497 


-6.0505 


-4.0594 


90 


1.4445 


-6.0499 


-6.0508 


-4.0592 



For J greater than 1/2 there are three open channels, which increases the 
magnitude of the numerical problem. In addition, there are more channels for 
a given value of jmax- To illustrate the convergence for higher values of J, we 
consider the J = 5/2"*" state for an incident laboratory energy of 2.0 MeV. For 
this case there are nine elements in the K-matrix. To show the symmetry of the 
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final K-matrix, we list the off-diagonal elements in a separate table. The indices 
1, 2 and 3 refer to the states with (Lq, Jo) equal to (2, 3/2), (2, 5/2) and (4, 7/2), 
respectively. Table ^ illustrates the convergence of the diagonal elements as the 
number of channels is increased, and Table |^ demonstrates the convergence for 
the off-diagonal elements. 

Table 6. Values of the diagonal elements of the K- 
matrix for J = f ^ p-d scattering at Eiah ~ 2.0 MeV 
evaluated using different numbers of basis states and 



pmax — 


90 fm. 








iVc 


J max 


Kii X 10^ 


K22 X lO'' 


K33 X 10=* 


82 


4 


-8.44217 


-1.37566 


-1.40145 


130 


6 


-8.45303 


-1.37703 


-1.40217 


178 


8 


-8.45504 


-1.37725 


-1.40284 


226 


10 


-8.45541 


-1.37725 


-1.40337 



Table 7. Values of the off-diagonal elements of the K-matrix for J ~ p-d scattering at 
Slab = 2.0 MeV evaluated using different numbers of basis states and pmax = 90 fm. 





J max 


K12 X 10^ 


7^21 X 10^ 


Kvi X 10-* 


K31 X lO"* 


K2S X 10^ 


K32 X 10'' 


82 


4 


3.41079 


3.40639 


-7.24681 


-7.40645 


-1.49972 


-1.54607 


130 


6 


3.40752 


3.40611 


-7.33861 


-7.40684 


-1.52924 


-1.54380 


178 


8 


3.40673 


3.40617 


-7.37622 


-7.40731 


-1.53775 


-1.54093 


226 


10 


3.40650 


3.40626 


-7.39359 


-7.40771 


-1.54067 


-1.53952 



Finally, in Tables ^ and |^ we show the convergence as the value of Pmax is 
increased. One can see that a value of 90 fm is sufficient for an accurate solution. 

Table 8. Values of the diagonal elements of 
the K-matrix for J — I"*" p-d scattering at 
Slab = 2.0 MeV evaluated using jmax = 10 and 
different values of pmax. 



Pmax 


Kii X 10=* 


K22 X 10" 


Kas X 10^ 


30 


-8.47082 


-1.37961 


-0.83198 


50 


-8.45433 


-1.37684 


-1.43016 


70 


-8.45570 


-1.37721 


-1.40318 


90 


-8.45541 


-1.37725 


-1.40337 



4 Conclusions 

A new numerical method for solving the large matrix equations in configuration- 
space three-body scattering has been shown to be an efficient procedure for 
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Table 9. Values of the off-diagonal elements of the K-matrix for J = f p-d scattering 

at Ei:i], = 2.0 M(A' oA'aluatt-d u.^iii^ Jwj,:^ = 10 and diffcrt-iit valuer of f>ui>ix- 



Pmax 


Ki2 X 10^ 


1^21 X 10^ 


Ki3 X 10* 


Ksi X 10" 


K23 X 10* 


K32 X 10* 


30 


3.41844 


3.41881 


-6.24969 


-7.56114 


-1.30255 


-1.56836 


50 


3.40606 


3.40580 


-7.43409 


-7.43557 


-1.55406 


-1.55135 


70 


3.40656 


3.40632 


-7.40409 


-7.41772 


-1.54418 


-1. -54286 


90 


:!. 10650 


:!.iO()2G 


-7.:!9:!59 


-7. 10771 


-l.r)10()7 


-1.5:!952 



obtaining accurate solutions for the realistic nucleon-deuteron scattering prob- 
lem. In addition, the convergence as a function of the number of channels and 
the matching radius has been demonstrated. 
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